# I load the needed libraries
library(dplyr)
library(scales)
library(GoFKernel)
library(mvtnorm)
library(gplots)
options(warn=-1)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
Loading required package: KernSmooth
KernSmooth 2.23 loaded
Copyright M. P. Wand 1997-2009
Attaching package: ‘gplots’
The following object is masked from ‘package:stats’:
lowess
I load the functions from the class file:
source("class_MCMC.R")
I define then the function that I want to use as output of the MCMCs:
# Function to sampled from: n-dim gaussian with chosen sigmas and centers
# posterior_g_inhom = function (theta) {
# sigmas = c(1:length(theta))
# centers = c(seq(length(theta), 1))
# product = 1
# for (i in 1:length(theta)) {
# product = product * exp(-(theta[i] - centers[i])**2/sigmas[i]**2)
# }
# return (product)
# }
cauchy2_gauss1 = function (theta) {
sigmas = 2.5
centers = 0.4
product = exp(-(theta[1] - centers)**2/sigmas**2)
product = product * (dcauchy(theta[2], -5, 2) + 4*dcauchy(theta[2], 8, 3))
product = product * (dcauchy(theta[3], -10, 2) + 4*dcauchy(theta[3], 10, 4))
return (product)
}
chosen_function = cauchy2_gauss1
Then I only have to determine the parameters for the initialization = the "hyperparameters" of the simulations
# The initial parameters are:
init = c(4, 4, 4)
std = diag(1, 3)
N = as.integer(1e5)
burn_in = as.integer(1e4)
print_step = as.integer(1e2)
# print_init = as.integer(1e3)
N_tot = N + burn_in
# For Haario:
epsilon = 0.001
# MVTNORM
# Evaluate then the MCMC
mcmc_g = random_steps_mvtnorm (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE)
# Selecting the sequence after the burn-in
mcmc_g = mcmc_g[burn_in:N, ]
# Plotting the results
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 77.91727 %
# MVTNORM GIBBS
mcmc_g = random_steps_mvtnorm_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 84.53848 %
# # SIMPLE ADAPTIVE
# mcmc_g = random_steps_simple (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
# gamma_function = gamma_series_exp, halved_step = burn_in)
# mcmc_g = mcmc_g[burn_in:N, ]
# show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
# # SIMPLE ADAPTIVE GIBBS
# mcmc_g = random_steps_simple_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
# gamma_function = gamma_series_exp, halved_step = burn_in)
# mcmc_g = mcmc_g[burn_in:N, ]
# show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
# HAARIO
mcmc_g = random_steps_haario (func_wanted = chosen_function, theta_init = init, n_samples = N_tot,
sigma = std, print_accept=TRUE, t_0 = burn_in, eps = epsilon)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 17.14091 %
Final mean = 0.3491998 4.911388 6.045531
Final covariance matrix =
[,1] [,2] [,3]
[1,] 6.292367 4.707983 12.64611
[2,] 4.707983 925.109283 102.26713
[3,] 12.646114 102.267133 3177.94219
# HAARIO GIBBS
mcmc_g = random_steps_haario_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot,
sigma = std, print_accept=TRUE, t_0 = burn_in, eps = epsilon)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 40.64273 %
Final mean = 0.4176557 4.912573 6.641108
Final covariance matrix =
[,1] [,2] [,3]
[1,] 7.149868 7.151321 10.16562
[2,] 7.151321 429.587025 127.00971
[3,] 10.165624 127.009711 721.32290
# RAO
mcmc_g = random_steps_AM_rao (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
gamma_function = gamma_series_exp, halved_step = burn_in/2)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 37.30636 %
Final mean = 0.2903034 5.115504 4.283763
Final covariance matrix =
[,1] [,2] [,3]
[1,] 3.1250070 -0.5717082 -0.1431036
[2,] -0.5717082 166.5288718 8.3458029
[3,] -0.1431036 8.3458029 758.3823423
# RAO GIBBS
mcmc_g = random_steps_AM_rao_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
gamma_function = gamma_series_exp, halved_step = burn_in/2)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 49.92303 %
Final mean = 0.4513963 8.014746 6.460089
Final covariance matrix =
[,1] [,2] [,3]
[1,] 3.49274003 -1.784855 0.03698277
[2,] -1.78485484 1174.422909 5.57069330
[3,] 0.03698277 5.570693 144.95974428
# GLOBAL
mcmc_g = random_steps_global (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
gamma_function = gamma_series_exp, halved_step = burn_in)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 45.98545 %
Final mean = 0.6499999 6.568418 5.424311
Final lambda = 0.7988679
Final covariance matrix =
[,1] [,2] [,3]
[1,] 3.351910 2.42679 -1.646531
[2,] 2.426790 222.62122 -12.976591
[3,] -1.646531 -12.97659 238.257206
# GLOBAL GIBBS
mcmc_g = random_steps_global_gibbs (func_wanted = chosen_function, theta_init = init, n_samples = N_tot, sigma = std, print_accept=TRUE, t_0 = burn_in,
gamma_function = gamma_series_exp, halved_step = burn_in)
mcmc_g = mcmc_g[burn_in:N, ]
show_results(mcmc = mcmc_g, init = init, std = std, step = print_step)
Acceptance rate = 73.8603 %
Final mean = 0.171423 9.719263 6.075724
Final lambda = -1.980399
Final covariance matrix =
[,1] [,2] [,3]
[1,] 3.7957458 -0.2153405 -1.980924
[2,] -0.2153405 329.5516120 6.426609
[3,] -1.9809240 6.4266093 334.293764
